On the relevance of numerical simulations to booming sand 
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We have performed a simulation study of 3D cohesionless granular flows down an inclined chute. 
We find that the oscillations observed in [L. E. Silbert, Phys. Rev. Lett, 94 098002 (2005)] near 
the angle of repose are harmonic vibrations of the lowest normal mode. Their frequencies depend 
on the contact stiffness as well as on the depth of the flow. Could these oscillations account for the 
phenomena of "booming sand"? We estimate an effective contact stiffness from the Hertz law, but 
this leads to frequencies several times higher than observed. However, the Hertz law also predicts 
interpenetrations of a few nanometers, indicating that the oscillations frequencies are governed by 
the surface stiffness, which can be much lower than the bulk one. This is in agreement with previous 
studies ascribing the ability to sing to the presence of a soft coating on the grain surface. 

PACS numbers: 45.70.Ht, 46.40.-f, 91.60.Lj 



A "booming" or "singing" dune is a sand dune that 
emits a loud sound when an avalanche occurs on its slip 
face. The sound can be very loud - audible up to 10 
km away - and has a well defined frequency of order one 
hundred Hertz. The physical origins of this phenomenon 
are still matter of debate in spite of much experimen- 
tal and theoretical work [lJ-Q - A successful theory of 
booming sand must explain why the sound is not a mix 
of a wide range of frequencies, and therefore several fre- 
quency selection mechanisms have been proposed. We 
now summarize them non-exhaustively. According to 
Andreotti for a given grain size, the frequency is 
set by the shear rate of the shear band separating the 
avalanche to the static part of the dune. This first ex- 
planation is not compatible with the one proposed by 
Douady et al. Q, in which the frequency is set by a 
resonance of the flowing layer. A third explanation has 
been proposed by Vriend et al. [||: The frequency is not 
selected by the properties of the avalanche but by the 
acoustical resonance induced by the stratification of the 
dune, explaining why the frequency may vary with time. 
More recently [2| , Andreotti and Bonneau show that the 
shear band separating an avalanche from the static part 
of the dune induces an amplification of guided clastic 
waves, leading to a linear instability. The frequency is 
then set by the maximum value of the instability growth 
rate. Finally Mills and Chevoir Q made the interesting 
remark that something like booming sand had already 
been observed in numerical simulations Q that exhib- 
ited spontaneous oscillations in granular flows near the 
angle of repose. These oscillations were interpreted as 
signs of intcrmittency near the jamming transition. 

Here, we examine Mills and Chcvoir's rc-intcrprctation 
of these oscillations as an acoustical phenomena - pos- 
sibly "booming sand" . We extend the simulations per- 
formed in Q by carrying out an extensive study of the 
influence of the contact normal stiffness and of the height 
of the flow. After a description of the method, we will 
show that the observed oscillations are multi-body har- 
monic oscillations. Their frequency is thus governed by 
the stiffness of the springs that model the repulsive inter- 



granular forces. We estimate a reasonable value for this 
quantity and discuss the relevance of the simulations. 

Method. We use our own implementation 0] of the 
classical "Discrete Element Method" method where New- 
ton's equations of motion for a system of N "soft" grains 
are integrated. This requires giving an explicit expression 
for the forces that act between grains. Such a technique 
is able to reproduce successfully experimental results for 
gravity driven flows 043, sheared systems [3], granu- 
lar materials close to jamming [l5[, silos [l6| or rotating 
drums fl7H20| . The Discrete Element Method is well 
known and can be found in many papers 0-[2l| . There- 
fore, we just present here the forces used in this work 
(and also in [7| ) . For the normal force between two over- 
lapping spheres we use a standard linear spring-dashpot 
interaction model [22|: f™ = k n S n — 7™v™, where S n is 
the normal overlap, k n is the spring constant, 7™ the 
damping coefficient and v™ the normal relative veloc- 
ity. The damping models the dissipation characteristic 
of granular materials. Likewise we model the tangential 
force as a linear elastic and linear dissipative force in the 
tangential direction: f* = — fc*5* — 7*v*, where k l is the 
tangential spring constant, 5* the tangential overlap, 7* 
the tangential damping and v* the tangential velocity at 
the contact point. The magnitude of 6 f is truncated as 
necessary to satisfy Coulomb law: |f*| < /i|f"|, where /i 
is the grain-grain friction coefficient. 

Numerical set-up. As in Q, we simulate gravity-driven 
chute flow. In every way the parameters are the same as 
in except we vary the normal stiffness between grains 
and the height of the flow. The grains are monosized 
(diameter d, mass m). Unless otherwise specified the 
number of grains is N = 8000. The chute consists in 
a 3D cell whose base is flat and rectangular with size 
20d x lOd It can be inclined relative to the horizon- 
tal by an angle 9 (angle between the horizontal and the 
long axis of the base) and is periodic in the directions 
tangent to the base. The bottom of the cell is obtained 
by pouring under gravity g a large number of grains in 
the cell (6 = 0°) and by fixing those that are in contact 
with the base. This disordered layer of fixed grains is 
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FIG. 1: (color online) (a) Average kinetic energy per particle 
for an angle of inclination 8 = 19.5° and k n d/mg = 2 x 10 6 . 
Fluctuations (main panel) at large time scales and oscillations 
at small time scales can be observed (inset), (b) Average 
kinetic energy per particle for 6 = 20° (the angle of repose is 
19.1°) and for different values of k n . 
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FIG. 2: (color online) (a) Frequency (cycles per time unit) of 
the kinetic energy oscillations (o), corresponding 1/2 power 
fit (dashed line), amplitude of the corresponding Fourier com- 
ponent (□) and corresponding (fc n ) _1 fit (full line). The angle 
of inclination is 6 — 20°. (b) same curves in a log-log plot. 



sufficient to prevent crystallization throughout the sys- 
tem. Before any measurements are taken, the inclination 
is increased to 8 w 30°, causing a rapid flow that erases 
any influence of the initial state. The angle 8 is then 
set to a final value that is cited in the captions or text 
below. The following values of the parameters are used: 
2 x 10 4 < k n d/mg < 2 x 10 6 , k* = 2k n /7, 7* = and 
/i = 0.5. The value of 7™ is adjusted to obtain a normal 
restitution coefficient e„ = 0.88 @. We use dimension- 
less quantities by measuring distances, times, and clastic 
constants respectively, in units of d, \fdfg and mg/d. 

Numerical Results. In accord with previous results Q , 
we observe an oscillation near the angle of repose which is 
defined, for a given flow height, by the angle below which 
the flow stops. This motion can be identified through 
measurements of the total kinetic energy. Note that in 
all our simulations the rotational kinetic energy is much 
lower than the translational one. As reported in Fig. [TJi, 
the kinetic energy displays regular oscillations (character- 
istic frequency ps \fgjd) and irregular fluctuations with 
lower frequency (ps 0.1 to Q.2y/g/d). We next turn our 
attention to the spring stiffness k n by studying its effect 
on the amplitude and on the frequency of the oscillations. 
As reported in Fig. [Tb, the oscillations indeed depend on 
the spring stiffness. As k n increases, their frequency in- 
creases while their amplitude decreases. Fig. [5] shows 
the frequency and amplitude of these oscillations versus 
spring stiffness. Frequency scales with vfc™, whereas am- 
plitude scales with 1/fc™. In the limit k n — > 00, the os- 
cillations disappear. However, the avalanche does not 
disappear in this limit, indicating that the avalanche and 
the oscillations are two separate processes. For this rea- 
son, we do not think that the oscillations should be in- 
terpreted as arising from the jamming transition, or as 
avalanche precursors. Now let us consider the oscillation 
frequency. It scales as 1 /t c where t c is the two-body col- 
lision time. Assuming that a two particle collision is a 
demi-cycle of a damped harmonic oscillator leads to t c = 



n/[2k n /m ~ (7 Tl ) 2 /m 2 ] 1/2 « it / \2k n /to) 1 / 2 cx l/Vfe". 
This scaling of the frequency suggests that the oscilla- 
tions are harmonic vibrations of a normal mode. Let us 
estimate the frequency of the vertical normal modes. In 
the simulations, there are a certain number nj, of lay- 
ers of grains resting on the bottom of the chute. We 
model the granular bed as a one dimensional chain of Hl 
masses connected by linear springs. This model is the lin- 
ear harmonic chain, used in solid state physics as a very 
elementary model of phonons (23[. The damping added 
to the particle interactions affects short wavelength vi- 
brations most strongly, but long wavelength ones only 
weakly [24j . Thus the motion is dominated by the longest 
possible wavelength. If the bottom of the chute is consid- 
ered to be fixed, and the top surface is free, the longest 
wavelength is four times the dep th of the lay er. This 
leads to a frequency f nL / y/g/d ~ (1 /4n£ )yVk n d/ mg. 
With riL = 40 and k n d/mg = 2x 10 6 , one obtains 
/40 ~ 9\/g/d. This value is twice as large as the one 
shown in Fig. [2] The difference between the prediction 
and the simulation is probably caused by the proximity 
of the stability threshold. The macroscopic stiffness of 
a granular packing decreases as a yield condition is ap- 
proached. Here, the yield condition (i.e the avalanche) is 
reached when 8 equals the angle of repose. To support 
this idea, simulations were carried out using exactly the 
same protocol but with a final angle of 6 = 0°. The sys- 
tem comes rapidly to rest, with the longest-lived motions 
being persistent oscillations at frequencies closer to the 
predicted values. But the central point is that the model 



reproduces the correct scaling with stiffness f nL ~ v^ n . 
Our results also predict that frequency should diminish 
as the granular layer is made deeper: f nL k 1/n^. To 
check our analysis we carried out numerical simulations 
with k n d/mg = 2 x 10 6 , 8 = 20° and for several number 
of grains (8000 < N < 32000). The results shown in 
Fig. [3] confirm the predicted scaling. We conclude, there- 
fore, that the oscillations observed in the simulations are 
simply harmonic oscillations of the lowest normal mode. 
Relevance of the numerical results. Up to this point, 
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FIG. 3: (color online) Frequency (cycles per time unit) of the 
kinetic energy oscillations (o) versus the height of the flow 
and corresponding fit (dashed line). The angle of inclination 
is 9 — 20° and the normalized stiffness k n d/mg = 2 x 10 6 . 



we have discussed our results solely within the framework 
of the very idealized model where a relatively small num- 
ber of perfect spheres flow down a fixed inclined plane, 
interacting via linear spring forces. The linear force law 
allowed us to reproduce the oscillations under exactly 
the same conditions as Silbert Q, and also facilitated 
the analysis of the observed frequencies. The idealized 
context of our work is emphasized by the exclusive use 
of dimensionlcss units. We did not try to give physical 
values to any parameter such as the particle diameter 
d. For the remainder of the paper, we will investigate 
the suggestion of Mills et al. |6( and discuss the rele- 
vance of our findings to "booming sand" . This means 
that we must assign physical values to all quantities. The 
most troublesome (and the most important) parameter is 
the spring stiffness k n . Real sand grains have non- linear 
force-displacement laws that cannot be characterized by 
a single spring constant. We will, nevertheless, try to 
straddle the difference between the model and physical 
system by choosing a single value of k n relevant to the 
oscillations. This approach may be problematic: Modi- 
fying k n affects the rheology of the flo w \)M , even for the 
highest values of k n used in this work [12[ • To determine 
the appropriate value of k n , we will consider a slightly 
less idealized model where the grains are spheres made 
out of an isotropic clastic material with Young modulus 
E and Poisson ratio v. The contact force between two 
such spheres is given by the Hertz law: 



EV2d 
3(1 -v 2 ) 



(1) 



We will take k n to be the stiffness seen by the grains when 
they oscillate at low amplitude about their equilibrium 
positions. We have thus k n = dF n /dS = (3/2)F n /S. 
Rearranging this equation, we obtain 

F n = 2k n 8/3. (2) 

In this equation, k n is not a constant: k n oc 5 1 / 2 . Wc 
seek a typical value of k n that will be determined by the 
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FIG. 4: (color online) The dimensionless stiffness of the con- 
tact as estimated by Eqs. ((4]) and ([5]) for a variety of diameters 
d (in meters) and number ni of layers of particles. 



typical contact force. We suppose that this force is of the 
order of the weight of a column of Ul grains: 



F n = rijjmg = nriLpgd 3 /6, 



(3) 



where p is the density of the material making up the 
grains. Eqs. ([I]), ©, and ([3]) involve the unknowns F n , 
k n , and 6. Combining them yields to 
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We will consider modifying Ul and d, while sweeping all 
the other constants into a single parameter that depends 
on the material: 
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Putting in values of appropriate for glass: g = 10m/s 2 , 
p = 2.4 x 10 3 kg/m 3 , v = 0.2, and E = 50 GPa, one 
obtains K = 2.3 x 10 4 m 2 / 3 . Various values of the di- 
mensionless stiffness k n d/mg, are shown in Fig. 2) Note 
that at d « 1 mm, even the softest contact (those that 
support the weight of only one grain) have a stiffness 
k n d/mg w 10 6 , i.e., equal to the half of the highest value 
in Fig.[TJ Diameters typical of "booming dunes" - around 
200 fj,m - lead to even stiffer contacts. Now let us see 
whether the oscillation frequencies in the simulation cor- 
respond to those of "singing sand" . The unit of frequency 
used here is \fgjd ~ 220 Hz for d = 200 pm (which is, 
as mentioned before, the typical diameter for grains of 
booming dunes). The observed frequencies in the field 
and in laboratory are around 90 Hz, i.e. 0Ay/g/d (see 
Fig. 3 of [3or table I of 0). According to laboratory ex- 
periments 3 , flow heights are of order of several centime- 
ters. For d = 200 pm, this corresponds to H/d w 500. 
Extrapolating the data in Fig[3]to H/ d = 500 leads to fre- 
quencies tantalizingly close to those observed. But Fig. [3] 
was obtained for a fixed stiffness k n d/mg = 2 x 10 6 , much 
lower than the value k n d/mg = 5.3 x 10 7 predicted by 
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Eq. ([5]) for Ul — 500 and d = 200 //m. If the oscillations 
are indeed the origin of booming sand, the effective stiff- 
ness of the grains must be much lower than predicted by 
Eq. ([5]). To explain this disagreement, let us determine 
the overlap 5 between two contacting grains from Hertz 
theory which assumes those grains as perfect spheres, 
without any surface asperities. Eqs. ([I) and ^ yield to 

^« 

For d — 200 /im and Ul = 500, we obtain 5 « 3 nm. This 
is a maximum value of S, concerning contacts that sup- 
port the weight of 500 grains. Contacts near the free sur- 
face of the flow (til < 5) have S < 2 A. For such overlaps, 
the contact between two real grains will be dominated by 
surface properties that might be quite different from the 
bulk properties considered in the Hertz model. In partic- 
ular asperities or a layer of silica gel (as proposed in (2f|) 
could significantly reduce the stiffness seen by acoustical 
waves. These speculations, however, can only be con- 
firmed by examining the sand grains themselves. These 
surface effects could also account for the rapid increase 
in sound speed with depth Q . The contact stiffness near 
the dune surface would be anomalously soft and fixed by 
grain surface properties, whereas deeper in the dune, at 
higher contact forces, the stiffness would be much higher 
and dominated by the bulk properties. 

Conclusion. We have examined the oscillations ob- 
served Q in numerical simulations of granular beds near 



the angle of repose. These oscillations arc harmonic 
vibrations of the lowest vertical mode of the bed, and 
their frequency obeys the expected dependency on par- 
ticle stiffness and bed depth. They are not part of the 
avalanche motion but may be connected to "booming 
dunes" , if the effective contact stiffness is about 20 times 
smaller than expected from the Hertz contact law. Such 
a reduction of stiffness is possible because the Hertz law 
predicts extremely small grain overlaps, indicating that 
the stiffness is dominated by surface properties instead of 
bulk ones. If the oscillations are indeed related to "boom- 
ing sand" , it would mean that the sound originates from 
a resonance inside a flowing layer, similar to the expla- 
nation presented in (3). However, several issues remain 
open. For example, what excites the oscillations? One 
clue is given by the effect of polydispersity: we performed 
simulations with uniform size distribution of width 2da, 
with a = 0.05, 0.1 or 0.2, and observe that the oscilla- 
tions disappear for a = 0.2, consistent with the obser- 
vation [261 ] that booming sand has a narrow grain size 
distribution. This suggests that the resonance is excited 
by quasi-periodic collisions in the shearing layer, consis- 
tent with most explanations that have been presented. 
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